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We study the field dependence of the antiferromagnetic spin-| Heisenberg model on the square lattice by 
means of exact diagonalizations. In a first part, we calculate the spin- wave velocity c, the spin-stiffness ps, 
and the magnetic susceptibility x± and thus determine the microscopic parameters of the low-energy long- 
wavelength description. In a second part, we present a comprehensive study of dynamical spin correlation 
functions for magnetic fields ranging from zero up to saturation. We find that at low fields, magnons are well 
defined in the whole Brillouin zone, but the dispersion is substantially modified by quantum fluctuations com- 
pared to the classical spectrum. At higher fields, decay channels open and magnons become unstable with 
respect to multi-magnon scattering. Our results directly apply to inelastic neutron scattering experiments. 

PACS numbers: 75.10.Jm 75.40.Mg 75.40.Gb 



I. INTRODUCTION 

The antiferromagnetic spin-^ Heisenberg model on the 
square lattice has been investigated in great detail over the past 
two decades, primarily because the parent compounds of su- 
perconducting copper-oxide materials are well described by 
this model On the square lattice, and more generally on 
any sufficiently connected bipartite lattice, the model is mag- 
netically ordered at zero temperature and the long- wavelength 
properties such as the magnetic order^, the low-energy excita- 
tions'^ or the temperature dependence of the magnetic correla- 
tions^'^'^ are well described within a semiclassical setting^'^'^^. 
At zero field, the staggered moment rris is however reduced 
by around 40% compared to its classical value, indicating the 
presence of sizable quantum fluctuations. One prominent ef- 
fect of fluctuations, the renormalization of the spin- wave dis- 
persion along the magnetic zone boundary between momenta 
(7r,0) and (7r/2,7r/2), was analyzed in Refs. 11,12,13,14,15. 
Their numerical and experimental findings showed that while 
the classical theory predicts a flat band along the boundary, the 
energy and the quasi-particle weight at (tt, 0) are substantially 
lower than at (7r/2 , 7r/2) . 

In the presence of a magnetic field, fluctuations are ex- 
pected to be even more important because of the additional 
interaction between spin waves and the field. Since the major- 
ity of studies of the Heisenberg model is devoted to the zero or 
small field situation, these effects are not yet well understood. 
The negligence of the strong field regime is not surprising - 
until recently, the prospects of experimentally probing the in- 
termediate or even high field physics were dire, given the huge 
exchange coupling of layered cuprate materials. This situa- 
tion has changed with the synthesis of CuBr4 and CuCU com- 
pounds^^'^^'^^ that have experimentally accessible saturation 
fields of less than 25 Tesla. In contrast, the high-Tc parent 
compounds have enormous saturation fields of several thou- 
sand Teslas. At low fields, spin- wave theory predicts only 
a weak renormalization of the magnon dispersion due to hy- 
bridization between single- and two-magnon states Zhito- 
mirsky and Chernyshev^^ showed that the effects of this hy- 
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FIG. 1: (Color online) Classical Heisenberg model in a magnetic 
field h directed along the z-axis. (a) For small fields, the spins are 
aligned antiferromagnetically in the xy -plane. The 0(3) symme- 
try of the model without magnetic field is reduced to 0{2) rotations 
within this plane, (b) For stronger fields, the spins develop a uniform 
component along the field and are thus canted out of plane. At the 
saturation field hs — 8JS, the spins are aligned ferromagnetic ally. 



bridization are more drastic at higher fields, leading to the in- 
stability of spin waves with respect to spontaneous decay. 

The aim of the present work is to provide a comprehensive 
numerical study of the Heisenberg antiferromagnet in mag- 
netic fields ranging from zero up to saturation. Rather than 
focussing on a particular feature, we aim to give an overview 
of the various properties induced by the coupling to the mag- 
netic field. Since our results for the dynamical spin structure 
factors are directly comparable to the inelastic neutron scat- 
tering cross section, we believe that this study is of interest to 
both theorists and experimentalists. 

The paper is organized as follows: In a first part, we dis- 
cuss static properties of the Heisenberg model in a mag- 
netic field, presenting results for the field dependence of the 
transverse susceptibility x± (Sect. Ill B), the spin-stiffness ps 
(Sect. IIIC), the spin-wave velocity c (Sect. HID), the static 
spin structure factors (Sect. HIE), and the staggered moment 
rris (Sect. IIIF). In Sect. IV, the main part of this work, we 
present our results for the dynamical spin correlation func- 
tions obtained by exact diagonalizations of clusters with up to 
64 sites. Our conclusion can be found in Sect. V. The reader 
primarily interested in the dynamical spin structure factors can 
skip the static part and focus on Sect. IV. 
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II. MODEL 

The isotropic spin- ^ Heisenberg model on the square lattice 
with nearest-neighbor interactions is described by the Hamil- 
tonian 



N=20, msO Y N=26, maO N=32, msO N=36, mal/6 N=40, ms1/4 



H 



(1) 



where is a spin-^ operator on site i, the summation (z, j) 
extends over nearest-neighbor pairs, and his 3. constant mag- 
netic field applied along the quantization axis. In what fol- 
lows, we set ^/iB 1, 5' being the gyromagnetic ratio and 
/i B the Bohr magneton. If not stated explicitly, energies are 
measured in units of the antiferromagnetic exchange coupling 
J > 0. For a system with N spins, the magnetization of a 
given state is defined as m = (Sf) /N. The fully polar- 
ized limit is thus reached at m = 1/2. 

In the absence of a magnetic field, the ground state of the 
classical model is the Neel state, in which all spins are aligned 
antiparallel with respect to one another, spontaneously break- 
ing the 0(3) symmetry of the Hamiltonian [Eq. (1)]. As the 
magnetic field is turned on, the symmetry of the model is re- 
duced to 0(2), describing rotations around the direction of 
the applied field. For small fields, the spins preferably align 
antiferromagnetically in the plane perpendicular to the mag- 
netic field, with a small uniform out-of-plane component, as 
depicted in Fig. 1. This uniform component becomes stronger 
with increasing field, until at the saturation field hg = 8JS 
all spins are aligned ferromagnetically. In the canted regime, 
the orientation of the spins can be decomposed into a stag- 
gered part perpendicular to the magnetic field and a uniform 
component directed along the field. 

The low-energy long-wavelength properties of the Heisen- 
berg model are well described by a nonlinear a model ^^'^^ 
whose Lagrangian density is defined as. 



C = 



(h-h X n)^ 



(2) 



Here n is a three-dimensional vector representing the orienta- 
tion of the staggered spin component subject to the constraint 
= 1, h is a constant magnetic field, and a dot denotes the 
time derivate. This model has three physical parameters: the 
spin-stiffness ps, the staggered moment rris and the uniform 
magnetic susceptibility x± in the direction perpendicular to n. 
The spin-wave velocity c is obtained from the hydrodynamic 
relation^^ 



Ps/X- 



(3) 



The a model description is valid at small fields, where the 
staggered moment is large. In addition to being renormal- 
ized by quantum fluctuations, see, e. g., Refs. 23,24, these mi- 
croscopic parameters also depend significantly on the strength 
of the magnetic field. The spin-wave velocity for instance de- 
creases with h and vanishes at the saturation field hs . We dis- 
cuss the field dependence of these three parameters in Sec. III. 

Finite-size effects in an 0(n) a model have been discussed 
in Refs. 25,26,27. For a quadratic sample with N sites, the 
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FIG. 2: ( Color online) Reciprocal lattices of the different finite -size 
clusters used in this work. Zero-field results are obtained from the 32- 
sites sample, while for the polarized regime, exact diagonalizations 
of clusters with upto N = 64 sites can be performed. The achievable 
magnetizations are indicated for every cluster Saturation is reached 
atm= 1/2. 



ground state energy density scales as 



eo(^') = eo - 1.437745 



2 J Ar3/2 
(n-l)(n-2) 
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PsN^ 



O 



, (4) 



allowing one to extract the spin- wave velocity c. Note that for 
n = 2, there is no contribution. The staggered magne- 

tization behaves as^"^'^^ 



ml{N) 



1 + 0.62075 



TV 



O 



(5) 

We will use these extrapolation formulas in Sec. Ill to estimate 
the spin- wave velocity and the staggered magnetization in the 
thermodynamic limit, taking into account the reduction of the 
symmetry from n = 3 to 2 in the presence of a magnetic field. 



III. STATIC PROPERTIES 

In this Sect., we discuss static properties of the Heisen- 
berg model. Readers more interested in the dynamical spin 
structure factors can directly proceed to Sect. IV. Our results 
are obtained by means of exact diagonalizations and quan- 
tum Monte-Carlo simulations using the SSE application^^'^^ 
of the ALPS library^^. The reciprocal lattices of the differ- 
ent clusters used in this work are shown in Fig. 2, together 
with the accessible magnetizations. Throughout this paper, 
we use the same color code for all finite- size data. Taking 
into account translation symmetries, the Hilbert space for a 
cluster of N sites at magnetization m encompasses roughly 
(Ar/2-mAr) states. Computations at low fields are most 
demanding, because of the huge Hilbert spaces requiring a 
significant amount of main memory. At zero field, exact diag- 
onalizations in the present study are limited to 32 or 36 sites, 
while in the high-field regime, we can include data from sys- 
tems with up to 64 sites. The largest Hilbert spaces involved 
encompass up to several 100 million states. Because the mag- 
netization m is a conserved quantum number, it is compu- 
tationally advantageous to work with fixed m rather than at 
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FIG. 3: (Color online) Inset: Uniform magnetization m as a func- 
tion of the applied magnetic field h. The dashed line is the classical 
result. Main panel: Difference to the classical curve Am(h) — 
m{h) — h/8. The solid black line represents the linear spin-wave 
approximation. Small red dots are obtained from the derivative of 
the extrapolated ground state energy, open squares represent ex- 
act diagonalization results from clusters with at least 40 sites, and 
open triangles indicate quantum Monte-Carlo results at tempera- 
tures T < 0.02 J. 



a given magnetic field h, as it is the case in experiments or 
in spin- wave calculations. To avoid inexact transformations 
between these two conjugate variables, we have chosen to 
present all our numerical results as a function of the magne- 
tization. Comparisons with spin- wave calculations are estab- 
lished by mapping the magnetic field h onto m via the inverse 
magnetization curve h{m) obtained within linear spin- wave 
theory. This choice is motivated by the fact that we prefer to 
focus on exact numerical results, keeping spin- wave approx- 
imations at the simplest level that qualitatively captures the 
physical properties. 



A. Uniform magnetization 



FIG. 4: ( Color online) Transverse susceptibility x± cls a function of 
the magnetization m and the applied magnetic field h. The classi- 
cal result (dashed line) is constant, while the linear spin-wave curve 
(solid line) obtained as the numerical derivative of the uniform mag- 
netization diverges in the limit m ^ 1/2 because of the logarithmic 
singularity in m{h). 

40 sites^^ and also quantum Monte-Carlo simulations^ We 
have reproduced these results in Fig. 3, together with our val- 
ues obtained from the derivation of the energy. The actual 
magnetization curve is presented in the inset, while the main 
panel shows the difference to the linear classical behavior. The 
overall agreement between linear spin- wave approximation^^ 
and numerical results is quite remarkable, the difference be- 
ing everywhere smaller than 5%. Our approach of using this 
magnetization curve to compare spin-wave results with nu- 
merical simulations is thus justified. The magnetization curve 
obtained from the extrapolated ground state energy is in ex- 
cellent agreement with Monte-Carlo simulations, especially at 
higher fields, where larger samples can be used for the extrap- 
olation. Close to saturation, the magnetization is expected to 
behave linearly with logarithmic corrections^^'^^'^^'^^, in con- 
trast to the root singularities encountered in one-dimensional 
systems^^'^^. 



On finite-size samples, the magnetization process occurs 
in finite steps. To obtain a continuous magnetization curve, 
one can either interpolate these steps along h and m or in- 
vert the derivative of the extrapolated finite- size ground state 
energy [Eq. (4)] with respect to the magnetization de^jdm. 
Using the finite- size scaling has the disadvantage that the en- 
ergies obtained from exact diagonalizations first have to be 
interpolated with respect to the magnetization m before they 
can be extrapolated to the thermodynamic limit because the 
magnetization steps of different samples only coincide for a 
few particular values of m. The interpolation is unproblem- 
atic, because the energy is a smooth function of the magne- 
tization. Subsequently, the extrapolation is performed for all 
values of m that are realized in at least one sample. The mag- 
netization curve of the Heisenberg model has been calculated 
before by exact diagonalizations of clusters containing at least 



B. Transverse susceptibility 

The magnetic susceptibility measures the variation of the 
magnetization with respect to a small external field applied 
in addition to the arbitrarily large field h introduced in Eq. (1). 
We can thus construct a susceptibility tensor as 

9{0|5"|0) 

Defining the terms parallel and perpendicular with respect to 
the direction of the staggered moment, the transverse suscepti- 
bility x± measures the variation of the uniform magnetization 
with respect to the applied field 

dm 
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FIG. 5: Spin-stiffness ps as a function of the magnetization m and the 
applied magnetic field h. The results have been obtained by quantum 
Monte-Carlo simulations^^'^^ using theALPS^^ application and spin- 
wave calculations. The initial increase of the spin- stiffness can be 
attributed to a reduction of quantum fluctuations, while the linear 
decrease at high fields is due to the canting of the spins towards the 
direction of the applied field. 



In Fig. 4, we plot the transverse susceptibility as a function 
of the magnetization and the magnetic field. The upper x-axes 
shown in Figs. 4-6 and 9 represent the linear spin- wave map- 
ping between the magnetization m and the magnetic field as 
a function of which spin- wave results are unambiguous. Note 
that the classical results are plotted using the classical mag- 
netization curve. The zero-field quantum Monte-Carlo result 
is taken from Ref . 24, which we have arbitrarily chosen as the 
representative of a large body of works encompassing many 
different methods. The spin- wave result has been derived 
in Ref. 38. Exact diagonalization and quantum Monte-Carlo 
data^^ at finite field are obtained by numerically differentiat- 
ing the magnetization curves shown in Fig. 3. In the classical 
limit, the susceptibility is constant, x j_ = 1/8, and the hydro- 
dynamic relation [Eq. (3)] is exactly satisfied for all magne- 
tizations. Quantum fluctuations modify the shape of x± con- 
siderably and lead to the divergence in the limit m 1/2, 
stemming from the logarithmic singularity in the magnetiza- 
tion curve. We note that the linear spin- wave approximation 
captures these effects extremely well and is in perfect agree- 
ment with Monte-Carlo data. 



C. Spin-stiffness 

The elastic energy required to twist a spin arrangement is 
proportional to the spin-stiffness p^, see Eq. (2). For a planar 
order parameter of the form 

n(r) = [cos(/9(r), sin(/p(r), 0] , 

that slowly rotates with a constant twist dip = Lp{Yi) — (p{rj) 
between adjacent sites and vj, the difference between the 
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FIG. 6: (Color online) Spin-wave velocity c as a function of the 
magnetization m and the applied field h. Symbols correspond to 
the slope of the dispersion extracted from finite- size clusters of size 
N. Small gray dots represent the thermodynamic values obtained by 
extrapolations according to Eq. (4). The dashed and the solid lines 
are spin-wave results. 

energy density of the twisted and the collinear spin arrange- 
ment is given by 

6e = eo{S(p) - eo{5(p = 0) ^ J ps {S(pf , 

where the twist is applied along both space directions. In 
principle, it is possible to calculate ps in exact diagonaliza- 
tions^^, but it is far easier and more accurate to obtain the 
spin-stiffness from quantum Monte-Carlo simulations as the 
global winding number fluctuations'^^. In Fig. 5, we plot ps as 
a function of the magnetization and the magnetic field. Previ- 
ous results at zero field are taken from Ref. 24. Given the sim- 
plicity of the linear spin-wave approximation, the excellent 
agreement with Monte-Carlo results is remarkable. Because 
the spin-stiffness is proportional to the elastic energy required 
to deform the collinear Neel order, it provides a measure of 
the ordering tendencies, or inversely, the quantum fluctuations 
present in the antiferromagnet. The initial increase of ps with 
the magnetization can be attributed to a reduction of quantum 
fluctuations. The almost linear decrease at high fields how- 
ever is due to the canting of the spins and has nothing to do 
with fluctuations because the exact result is almost identical 
to the classical curve. The reduction of quantum fluctuations 
at low fields can also be seen in the staggered magnetization 
discussed in Sect. Ill F. 



D. Spin-wave velocity 

There are two ways to extract the spin- wave velocity c from 
exact diagonalization results: One can either determine the 
slope of the dispersion in the vicinity of the ground state mo- 
mentum, or, one can fit the finite- size scaling of the ground 
state energy to Eq. (4) and determine c from the fitting param- 
eters. This latter method has the disadvantage that the energies 
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FIG. 7: (Color online) (a) Finite-size extrapolations (dashed lines) 
of the energy eo according to Eq. (4) for magnetizations m = 0, 
1/8, and 1/4. The spin-wave velocity c is equal to the slope of the 
extrapolations. The error bars shown in Fig. 6 correspond to the 95% 
confidence intervals of the above fits, indicated by the thin solid lines, 
(b) Finite-size extrapolations of the staggered moment rnl according 
to Eq. (5). The error bars represent the 95% confidence levels of the 
linear regressions, also shown in Fig. 9. Zero-field data from clusters 
with more than 32 sites are taken from Ref. 31. 



first have to be interpolated with respect to the magnetization 
m, before the extrapolation can be performed, as discussed 
in Sect. Ill A. In Fig. 6, we show the results from both ap- 
proaches and compare them to spin- wave calculations^^. 

While the three methods give qualitatively similar results, 
the quantitative differences are non-negligible, especially in 
the low-field regime. Surprisingly, the zero-field spin-wave 
velocity extracted from the sample with 32 sites is in per- 
fect agreement with quantum Monte-Carlo data^"^. The re- 
sults from the 20 and 26-sites cluster are somewhat lower, 
suggesting a thermodynamic value slightly above the one ob- 
tained from the 32-sites system. The abrupt drop of the veloc- 
ity upon turning on the magnetic field is unphysical and the 
pronounced finite- size effects observed in this regime indicate 
that the actual value of c is substantially larger than the 32- 
sites cluster result. For increasing field, these effects become 
smaller and one could conclude that at m = 0.2, the thermo- 
dynamic value should be very close to the 36-sites result. This 
is however not quite correct because for any system size, the 
distance to the gapless mode is at least 7r/4, see Fig. 2 for the 
location of the momenta in different clusters, and thus still far 
from the cone center. For small and intermediate fields, one 
therefore systematically underestimates the spin-wave veloc- 
ity. At very high fields, the direction of the finite- size effects 
is reversed because the dispersion starts to resemble the ferro- 
magnetic eic cx k'^ form and the finite-size momenta are too far 
away from the gapless mode to capture this behavior. This is 
best seen in the longitudinal structure factors shown in Fig. 1 1 . 

Given the pronounced finite- size effects at small fields, one 
can argue that an extrapolation of the ground state energies ac- 
cording to Eq. (4) would give a more accurate value of c. On 
the contrary, at zero field, this method clearly underestimates 
the spin- wave velocity, as already noted in Ref. 41, while for 
small m, it largely overestimates c. We think that the huge ini- 
tial increase is due to the change of symmetry: While the ex- 
trapolation formula is modified when turning on the magnetic 
field, from n = 3 to 2, see Eq. (4), the finite-size clusters are 
too small to reflect this change. This interpretation is consis- 



tent with the observed overestimation of almost a factor 3/2. 
A reason for the underestimation of the zero-field spin-wave 
velocity might be the influence of higher order terms in the fit- 
ting formula [Eq. (4)], which forn = 3 also has a contribution 
1/7V^. Including higher order corrections leads to an increase 
of the velocity c, but also to inconclusively large error bars. 
We believe that the bump observed between m = 0.2 and 
0.3 is an artifact of the finite- size extrapolations, as it seems 
driven by the 36 and 40- sites cluster. To give an idea of the 
variability of the extrapolations, the fits for m = 0, 1/8, and 
1/4 are shown in Fig. 7(a). The 95% confidence intervals are 
shown as thin solid lines. They are also indicated in Fig. 6. 
Concerning more advanced spin-wave calculations of the dis- 
persion, we would like to refer the reader to a recent work by 
Kreisel et al^^, discussing a non-analytic behavior of the spin- 
wave velocity induced by the magnetic field. The validity of 
the hydrodynamic relation in a magnetic field has been con- 
firmed very recently in Ref. 38 within a spin- wave approach. 
At zero field, Hamer et al^^ verified this property using series 
expansions and spin- wave approximations. 



E. Static spin structure factor 

The static spin structure factor 



with a, /3 = X, +, — , is the Fourier transformation of the 
equal-time spin- spin correlators and can be directly measured 
in elastic neutron scattering. Fig. 8 shows the longitudinal 
and transverse structure factors obtained from exact diagonal- 
izations, see also Ref. 44, and spin- wave calculations^^ for 
several magnetizations along highly symmetric points in the 
Brillouin zone. Because there is no spontaneous symmetry 
breaking in finite systems, the exact diagonalization results 
for S^^{q) and 5'^^(q) are identical. Hence, it is useful to 
define a transverse component S^^ as 

s^y{k) = \[s+-{k) + s-+{k)] . 

In contrast, the spin-wave results are derived in the symme- 
try broken phase where the staggered spin component is di- 
rected along the x-axis. The transverse correlations shown in 
Fig. 8(b) represent the average 25^^ (k) = 6'^^(k) + S'^^(k). 
At zero field, the structure factors reflect the antiferromagnetic 
order with a pronounced peak at X = (tt, tt). As the mag- 
netic field is turned on, this contribution rapidly decreases in 
the longitudinal channel but remains dominant in the trans- 
verse structure factor, indicating staggered order in the plane 
perpendicular to the magnetic field. With increasing field, 
the spins are gradually canted out of the xy-plane and de- 
velop a uniform component that manifests itself at momentum 
r = (0, 0) in the longitudinal structure factor. At high fields, 
both curves are almost flat, except at the before mentioned 
two points. We note that linear spin- wave results capture the 
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FIG. 8: (Color online) Static longitudinal (a) and transverse (b) spin structure factors S^^{q) for various magnetizations m along a path 
connecting highly symmetric points in the Brillouin zone. Symbols represent exact diagonalization results, with numbers indicating the system 
sizes. Larger symbols stand for bigger clusters. Solid lines show the spin-wave results derived in Ref. 20. 



essential features qualitatively correctly, apart from the uni- 
form contribution to the longitudinal structure factor. Our re- 
sults are in good agreement with recent quantum Monte-Carlo 
data45. 



F. Staggered magnetization 

Similar to the uniform magnetization, the staggered mo- 
ment rris can be obtained as the derivative of the ground state 
energy density with respect to the staggered field. The linear 
spin- wave result has been derived in Ref. 46. Alternatively, 
rua can be obtained from the transverse structure factor as 



1 

TV- 




2 (/i = 0) 
1 (/i > 0) 



where N' is a normalization factor. On finite systems of size 
the static structure factors satisfy the relation 

^+-(k) =^-+(k) + 2m, 

contributing a term 1/N' io ml. This additional term makes 
the extrapolation according to Eq. (5) more difficult because 
one needs to take into account an additional parameter. To 
circumvent this problem, we only consider the S ^(k) con- 
tribution to the transverse structure factor. In the thermody- 
namic limit, both choices of the normalization N' are iden- 
tical, but on the small samples we study, there are sizable 
differences. The zero-field choice was introduced in Ref. 47 
and successfully applied to the J1-J2 model on the square 
lattice"^ ^. In the presence of a magnetic field, using only the 
S ^ (k) contribution which vanishes at k = F naturally leads 
to a normalization N' = N — 1. In Fig. 9, we show the 
staggered magnetization obtained from extrapolations based 
on Eq. (5), quantum Monte-Carlo simulations and spin-wave 
calculations. At zero field, quantum fluctuations are strong, 
leading to a considerable reduction of the ordered moment. 
In this case, all methods give a similar value slightly above 



mg ^ 0.3. As the field is turned on, rris increases and the fluc- 
tuations thus become less important, giving rise to an almost 
classical regime, in which the dispersion is well described by 
linear spin- wave theory. This aspect will be discussed in more 
detail in Sect. IV. After the maximum dXm ^ 0.15, the stag- 
gered moment is again reduced because of the canting of the 
spins and eventually vanishes in the saturation limit. While 
the extrapolation of exact diagonalization data according to 
Eq. (5) gives a very good result at zero field^^''^^ it consider- 
ably underestimates the staggered moment at small fields. To 
get an idea of the variability of the extrapolations, fits at mag- 
netizations m = 0, 1/8, and 1/4 are shown in Fig. 7. The 
95%-confidence intervals are shown as thin solid lines. They 
are also indicated in Fig. 9. These confidence intervals should 
not be interpreted in a strict statistical sense, but are rather in- 
tended to remind the reader of the fact that the extrapolations 
to the thermodynamic limit are based on only a handful of 
small clusters. The inaccuracy at small fields is related to the 
smallness of the clusters used for the extrapolations as well as 
the lack of analytical predictions for higher order finite- size 
corrections. Similar to the extraction of the spin- wave veloc- 
ity, the method is much more accurate at higher fields, where 
bigger samples can be taken into account. In contrast to finite- 
size extrapolations, spin- wave results are in excellent agree- 
ment with quantum Monte-Carlo data obtained on a cluster 
with A/" = 32 X 32 sites. 



IV. DYNAMICAL PROPERTIES 

In the main part of this paper, we study the dynamical spin 
correlation functions 

1 I \ 

id ■'-'^ 

(6) 

where, as before, N is the number of sites on the lattice. This 
quantity is directly measurable in inelastic neutron scattering 
experiments. In what follows, we consider a (3 = zz, which 
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FIG. 9: (Color online) Staggered magnetization rris obtained from 
finite-size extrapolation according to Eq. (5), quantum Monte-Carlo 
simulations, and spin-wave calculations. 

we refer to as longitudinal correlations because they are par- 
allel to the applied field, and af3 = H — , — h, allowing us to 
construct the transverse correlation functions defined as 

S-y{LO,q) = ^[S+-{lo,ci) + S-+{LO,q)] . 

In exact diagonalization one can rather easily calculate 
spectral functions directly in real frequency, using the so 
called continued-fraction technique^^'^^. Combined with the 
Lehmann representation 

S''^{iO,q) = ^ ^''(q) \^^)\^ 5{L0-[Ei-E^]) , (7) 

where ^rn is the ground state at a given magnetization m with 
energy Em and are excited states with energies Ei, this ap- 
proach allows one to obtain unbiased exact results, including 
residues, although only on finite- size samples. 

The Zeeman term —h ^ • 6'f = —hS^^^ is usually not in- 
cluded in exact diagonalization calculations because it is con- 
stant within a given S^^^ sector. While this omission does 
not affect the longitudinal structure factors, extracting the fre- 
quency dependence of the transverse spin correlations is more 
subtle because the two contributions act between different S^^^ 
sectors. This is best understood by tracing back the origin 
of the poles in the dynamical structure factors to the eigen- 
levels in the tower of states. In Fig. 10(a), we present part of 
the tower of states for the 3 2- sites cluster for 6 < < 10. 
The top panel shows the spectrum of the Hamiltonian without 
magnetic field, while in the bottom panel, we have subtracted 
the contribution hS^ , for a magnetic field h corresponding to 
m = 1/4, i.e., 5*^ = 8 on this cluster. Realigning the tower in 
the presence of a magnetic field is thus equivalent to shifting 
the spin lowering and raising contributions to the transverse 
structure factors in such a way that the lowest-lying poles of 
the k = Q modes coincide. The most important properties 
of the tower of states are the alternation of the ground state 



momenta between (0, 0) and (tt, tt) in adjacent sectors and 
the proportionality of the ground state energies io S^{S^ ^1). 

Figs. 10(b) and (d) show the longitudinal and transverse 
structure factors obtained from this spectrum. In Fig. 10(c), 
we have zoomed in on a part of the tower of states to facilite 
identification of eigenlevels and poles. As expected, at a mag- 
netization m = /N , the main magnon branch of the longi- 
tudinal correlations is made of spin S = states. Similarly, 
in the transverse channel, the magnon dispersion originates 
from spin S = ^1 states. As an exception to this rule, the 
pole at X = (tt, tt) in the longitudinal structure factor origi- 
nates from the lowest-lying spin S = ^1 level, as can be 
seen in Fig. 10(a). Its energy is equal to the magnetic field. In 
this sequence of pictures, it is easy to see that while the low- 
est lying poles of the spin raising and lowering contributions 
to the transverse structure factor coincide, higher lying levels 
are slightly shifted with respect to one another. 



A. Synthetic view of the dynamical structure factors 

In exact diagonalizations, the restriction to relatively small 
systems makes it difficult to extract wave vector dependent 
quantities to the thermodynamic limit. Only highly symmet- 
ric points like M = (0,7r) or = (7r/2,7r/2) are common 
among several systems sizes up to 64 sites. By construction, 
r = (0,0) and X = (tt, tt) are present in any cluster, but the 
physics at these points is usually quite simple. In the presence 
of a magnetic field, bigger clusters can be used for the calcu- 
lations, but the main difficulty lies in the fact that the magne- 
tization steps are in general not compatible because they also 
depend on the system size. It is therefore basically impossi- 
ble to present an exact superposition of all dynamical correla- 
tion functions obtained from different systems sizes at a given 
magnetization. However, since the shape of the dynamical 
structure factors does not change substantially from one mag- 
netization step to another, it is nevertheless possible to con- 
struct quite accurate interpolations. For a given magnetization 
m and a cluster of size A^, we chose to either take the ex- 
act data if available or instead shift the whole structure factor 
obtained for m' closest to m to the position determined by lin- 
ear interpolation between the lowest-lying poles at ± l/A'. 
While this interpolation procedure does not yield exact results, 
it allows one to visualize data from different clusters in a sin- 
gle plot, thus providing a nice overview of the observed fea- 
tures. In Fig. 11, we show such a synthetic superposition of 
the longitudinal structure factors obtained for various magne- 
tizations along a path through the Brillouin zone. The zero- 
energy mode of the longitudinal structure factor at momen- 
tum (0, 0) is proportional to the magnetization squared, i.e., 
S^^ico = 0, q = F) = (5'^)^. In what follows, we omit this 
contribution. Note that the magnetization steps in the first row 
are not equidistant. At low fields, m < 1/6, we exclusively 
use the 32-sites cluster for all calculations. 

At zero field, the ground state of the infinite systems spon- 
taneously breaks the 0(3) rotational symmetry of the Hamil- 
tonian [Eq. (1)], leading to the existence of two Goldstone 
modes. The distinct pole at k = X collapses to the ground 
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FIG. 10: (Color online) Tower of states (a,c) and longitudinal (b) and transverse (d) structure factors obtained from the 32-sites cluster at 
magnetization m — 1/4, corresponding to — 8. The ground state energy of the = 8 levels is set to zero. The location of the poles in 
the structure factors can be traced backed to the energy levels in the tower of states. To facilitate identification, the strongest signatures are 
indicated by the symbols of the corresponding eigenlevels. In the presence of a magnetic field, constructing the transverse structure factor is 
not straightforward because the relevant towers for and S ^ start at different energies. The alignment procedure is illustrated in the 
bottom panel of (a). 



state energy as because it originates from a spin-1 level 
that is part of the Anderson tower, see Refs. 8,47. The strong 
intensity reflects the dominantly antiferromagnetic character 
of the spin alignment. An interesting zero-field property dis- 
cussed in detail in the literature^ and mentioned in 
the introduction, is the dispersive feature observed along the 
magnetic zone boundary between M and S. In Fig. 12, we 
present a finite-size scaling (solid lines) of the spin-wave en- 
ergies and the quasi-particle residues at momenta (7r/2, 7r/2) 
[N = 16, 32] and (tt, 0) [TV = 16, 20, 32, 36]. The extrapola- 
tion to the thermodynamic limit is in surprisingly good agree- 
ment with quantum Monte-Carlo simulations^^ and series ex- 
pansion calculations^^. 

Let us now return to the discussion of Fig. 11. A finite 
magnetic field reduces the symmetry to 0(2), which is again 
spontaneously broken in the ground state, leading to only one 
gapless mode. The second mode now has a gap proportional to 
the magnetic field. The location of the distinct pole at k = X 
thus no longer scales to zero in the thermodynamic limit, but is 
simply equal the magnetic field. As mentioned in Sect. Ill F, 
discussing the transverse moment rris, quantum fluctuations 
initially decrease with increasing magnetic field and reach a 
minimum at around m ~ 0.15. It is therefore not surprising 
that the classical dispersion is almost identical to our exact 
results in the vicinity of these polarizations. Comparing the 
zone-boundary dispersions (from M to S) for different mag- 
netizations shown in the first row of Fig. 11, one finds that 
the zero-field effect is inverted and that for m > 1/8, the 
spin- wave energy at (7r,0) is higher than at (7r/2,7r/2). At 
m = 1/8, for instance, the energy differs by ~ 5%, with a 
slope opposite to the one at zero field. 

With increasing magnetization, the weights of the magnon 
poles decrease, while at the same time, poles at higher ener- 
gies become more important. This modification is first ob- 
served at momenta close to X = (tt, tt), but spreads quickly 
along almost the whole path shown in Fig. 11. Only momenta 
around the gapless point F = (0,0) have well defined spin- 



wave poles for all meaningful magnetizations. This is the 
finite- size manifestation of the magnon instability predicted in 
Ref. 20. In the decaying regime, the lowest-lying poles have 
very small weight and are thus difficult to identify. As a guide 
to the eye, we have connected them in the case of the 64-sites 
cluster, see last row of Fig. 11, showing the presence of multi- 
ple spin- wave continua below the primary spin- wave branch. 
The possibility to resolve poles with very little weight in the 
spectral functions is a conceptual advantage of exact diagonal- 
ization compared to other methods like quantum Monte-Carlo 
simulations. Using this information, it is possible to detect the 
availability of phase space for the decay of single spin waves, 
even in the case where the lifetimes are large. 

Since our calculations were done on finite clusters, one 
must be careful interpreting results at fields very close to sat- 
uration, which might be obtained from small Hilbert spaces 
with only one or two flipped spins. These situations are clearly 
not representative for the infinite system and we thus only con- 
sider data from sectors with magnetizations m < mmax = 
1/2 -4/N. 

The magnon dispersion at h = 3.5 is in very good agree- 
ment with recent quantum Monte-Carlo simulations combined 
with stochastic analytical continuation"^^, compare second plot 
in the last row of Fig. 11 with Fig. 4 of Ref. 45. Concern- 
ing the magnon decay, our exact diagonalizations indicate that 
spin waves close to (tt, tt) have a very long, but finite lifetime. 
In contrast, the density plot presented in Ref. 45 suggests that 
spin waves are stable in this region. Around (7r,0), the sit- 
uation is much clearer, and both methods predict a region of 
instability, despite the fact that the two-magnon decay kine- 
matics of the spin-wave theory do not allow for an instabil- 
ity at this momentum. It would be interesting to investigate 
these subtle differences in more detail and compare the pred- 
ications of these two numerical approaches to more advanced 
spin- wave calculations^^. 
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FIG. 11: Synthetic superposition of the longitudinal dynamical structure factors along a path of highly symmetric points in the Brillouin 
zone. Different colors represent data from different clusters and the area of the symbols is proportional to 5*^^(0;, q). Dashed lines show 
the dispersions obtained within linear (harmonic) spin-wave theory and the solid line represents spin-wave results with first-order corrections 
(only shown for m = Oj. For magnetizations around m ^ 0.15, quantum fluctuations are almost negligible and the spin-wave dispersion is 
in good agreement with numerical results. At higher flelds, m > 0.3, fluctuations are again important and lead to the spontaneous decay of 
magnons. This instability is reflected in a reduction of weight in main peak accompanied by the appearance of small poles at lower energies. 
This process starts around q = X and spreads over almost the whole Brillouin zone. 
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FIG. 12: Finite-size scaling of the spin-wave energies (a) and the 
quasi-p article residues (b) at momenta (7r/2, 7r/2) and (tt, 0). 



B. Quantum fluctuations and finite-size effects 

After this overview, let us now look at some of the discov- 
ered features in more detail, using the exact dynamical struc- 
ture factors without any interpolation presented in Fig. 13. A 
false-color density plot of the longitudinal structure factor is 
shown in Fig. 14 after artificial broadening by an imaginary 
part T] = 0.05. Because data from all available clusters is 
superimposed, some of the features abruptly appear as more 
Hilbert spaces become accessible for computations. In gen- 
eral, we expect the longitudinal spin correlations at momen- 
tum k to resemble the transverse ones at momentum k + Q, 
i.e., 

S^^iLO, q)/5^^(q) « S^y{u, q + Q)/^^J'(q + Q) . 

This follows from the tower of states shown in Fig. 10. Given 
the similar structure of levels in neighboring sectors, it reason- 
able to expect the dynamical structure factors to reflect these 
similarities. The results displayed in Fig. 13 nicely support 
this hypothesis. 

At zero field, it is well known that the dispersion obtained 
within lowest-order spin-wave theory is renormalized by 
quantum fluctuations, leading to energies roughly 16% higher 
than this simplest prediction. As explained in Sect. Ill F, dis- 
cussing the ordered transverse moment rus , quantum fluctua- 
tions are minimal at magnetizations around m ^i:^ 0.15. This is 
exactly the regime in which the lowest order spin- wave disper- 
sion is in good agreement with exact numerical calculations. 
However, away from this point, there are sizable deviations. 
Similar to the zero-field situation, the most pronounced effect 
is observed at the magnetic zone boundary along (tt, 0) and 
(7r/2, 7r/2), for which spin- wave calculations predict a flat en- 
ergy = 2 J, independent of the magnetization. In contrast, 
our field dependent results at (7r/2, 7r/2) [and to a lesser ex- 
tent also (tt, 0)] reveal a large negative slope, not at all com- 
patible with the flat linear spin- wave results. Away from the 
magnetic zone boundary, the qualitative features of the field 
dependence are nevertheless nicely captured by the harmonic 
spin-wave spectrum. Going beyond linear spin-wave theory, 
one can include the first 1/ S corrections in perturbation the- 
ory, see Eq. 9 in Ref. 20. In Fig. 13, these perturbative cal- 
culations are shown as solid lines. As explained in Ref. 20, 
such a perturbative approach breaks down at larger fields, at 



which one should take into account the renormalized magnon 
energies and solve the Dyson equation self-consistently. For 
this reason, the solid lines do not extend over the whole field 
range. Given the rather good agreement between our exact 
numerical results and these perturbative energies, we conclude 
that the dominant dispersive features are well captured includ- 
ing 1/ S corrections. 

Our findings corroborate recent results revealing that non- 
collinear magnetic structures encounter sizable renormaliza- 
tion of the excitation spectrum when improving the linear 
(harmonic) spin- wave theory by the first 1/ S corrections^ ^ 
While this happens already in zero field for the triangular lat- 
^1^^52,53 ^j^j^ noncollinear 120° magnetic order, on the square 
lattice^^, this effect is driven by the magnetic field, and can 
thus be tuned in experiments. 

In the longitudinal structure factors away from (0,0) and 
(tt, tt), one observes a sometimes sizable jump of the energy 
of the first pole as the system is polarized. This discontinuity 
is caused by different finite- size effects present at zero field 
and for m > 0. The longitudinal magnon mode at 5^ = 
is in fact the 5'^=0-component of a spin 1 level. At m = 0, 
the energy thus contains a finite- size contribution of the An- 
derson tower that scales away as 1/A^. In contrast, at m > 0, 
the longitudinal structure factor couples predominantly states 
of equal spin, and the leading 1/A^ finite-size corrections are 
absent. 

Apart from the main spin-wave branch, a second distinct 
feature of the dynamical structure factors is the almost linearly 
increasing mode at low fields. Its intensity fades away as the 
field is increased and it is visible at all momenta. This is not 
a new mode, but a finite-size hybridization between the lon- 
gitudinal and the transverse spin structure factors, which can 
be explained by the Wigner-Eckart theorem. At a magnetic 
field corresponding to a magnetization m, the ground state 
is in a total spin Sm = mN state with polarization = Sm- 
The scattering operator entering S ^ couples only to to- 
tal spin S = Sm + 1 states in the spectrum because it is the 
+1 component of a spin 1 tensor operator. The scatter- 
ing operator entering S^^ couples to total spin S = Sm and 
S = Sm + 1 States in the spectrum, while S~ couples to total 
spin S = Sm — 1, Sm and Sm + 1 states. Using the Wigner- 
Eckart theorem one can prove that the part of the 5*^ -operator 
that couples to l^*^ + 1, Sm) states is actually related to the 
action of 5'+, and the following relation holds 



S''{ci,Lu)\s=s^+i=C-+ X S-+{q,LU + h), 



with 

c- 



+ _ 



(1,0; Sm 7 Sm I Sm + 1 7 Sm ) 



(I7 1 5 ^mj Sm\Sm + I7 Sm + 1) 



1 



a squared fraction of Clebsch-Gordan coefficients. As one in- 
creases Sm = Nm, either by increasing m at fixed system 
size N, or by increasing N at fixed m, the proportionality 
factor decreases as l/Sm for large Sm, and therefore, 
the shifted shadow of 5* ^ fades away, unless S ^ is suffi- 
ciently divergent. Similarly, one can show that S~^~ contains 
a shadow of S^^ which is generically Clebsch-Gordan sup- 
pressed at large m and N. 
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(b) Transverse dynamical structure factors 1/2 S^"' '""'^((ja,k) 
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FIG. 13: Longitudinal (a) and transverse (b) dynamical structure factors at selected momenta. Different colors correspond to different 
clusters and the area of the symbols is proportional to the correlation functions. Dashed lines represent the dispersion obtained within 
linear (harmonic) spin-wave theory, while solid lines are spin-wave results with perturbative 1/S corrections. The upper x-axes indicate the 
magnetic field obtained within linear spin-wave theory. Crosses represent zero-field series expansion re suit s^^ . A false-color density plot of 
the longitudinal spin correlations is shown in Fig. 14. 
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FIG. 14: (Color online) Density plots of the longitudinal dynamical structure factors S^^(k,co) for various momenta k. These plots are 
obtained from the raw data presented in Fig. 13(a) after artificial broadening by an imaginary part rj — 0.05. Structure factors from all 
available clusters are superimposed. Dashed lines represent the dispersion obtained within linear spin-wave theory and the upper x-axes 
indicate the mapping of the magnetization onto the magnetic field. 



C. Magnon instabilities 



We have seen that the location of the dominant pole in the 
longitudinal structure factor at momentum (tt, tt), and hence 
at (0,0) in the transverse spin correlations, is equal to the 
magnetic field. The corresponding plots in Fig. 13 thus rep- 
resent nothing but the inverse magnetization curve m{h). In 
the longitudinal structure factor, the lowest-lying pole at small 
magnetizations is also the dominant one, but interestingly, for 
m > 9/32, we observe an increasing number of poles with al- 
most negligible weight appearing below the main spin-wave 
feature. Because of their low intensity, we have indicated 
the locations of these poles by uniformly sized open sym- 
bols in the area denoted "zoom". Only weights larger than 
10~^ have been included in the plot, but most symbols repre- 
sent poles with weights around 10~^. For small clusters, i.e., 
N < 36, this group of levels first appears around m ^ 0.1, 
still above the spin-wave band, remaining roughly constant 
at energy uj ^ 3, until they cross the main magnon branch 
and lose some of their intensity. Data from bigger clusters 
suggests that these low-lying excitations form a continuum. 
Hence magnons in the main branch seem to have a very long 
but finite lifetime. Given the tiny weight of the low-energy 
excitations, one can conclude that while there is phase space 
available at these energies, the decay matrix element is ex- 
tremely small. The dominant pole in the longitudinal (tt, tt) 



structure factor in sector = Sm originates from the local 
copy of the ground state level in the neighboring spin sector 
with spin S = Sm + 1- In contrast, the almost negligible 
poles can be traced back to low-lying spin S = Sm states 
whose energies decreases with increasing magnetic field. The 
transverse structure factor at (0,0) does not exhibit such a 
continuum of low-intensity poles below the spin- wave branch 
because the S~ operator acts as a global spin lowering opera- 
tor creating an exact excited eigenstate, which gives rise to a 
single delta function in S'^^(q = F, cj). 

Within linear spin-wave approximation, magnons are sta- 
ble at any momentum and magnetic field, simply because the 
Hamiltonian does not contain terms that would allow a single 
magnon to decay into two spin waves. Knowing that such pro- 
cesses are however present in the original Hamiltonian, one 
can use purely kinematical arguments to study the necessary 
conditions for magnon instabilities. Following Ref. 20, we 
expand the spin-wave dispersion in the vicinity of the gapless 
mode at (tt, tt) leading to 

eq ^ cq[l^ aq q'^] , 

where |q — Q| <C 1 and c is the spin- wave velocity. The 
curvature aq is wave-vector dependent. Along the diagonal, 
one finds a = — 7/3 + 2/\^l — h'^, which changes sign at a 
critical field he = 2/V7hs ^ 0.755 hs, corresponding to a 
magnetization rric ~ 0.33. This signals the onset of the spin- 
wave instability^^. 
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(a) Interpolated weight of the first pole in the longitudinal dynamical structure factor 
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FIG. 15: Interpolated weights of the first poles in the longitudinal (a) and transverse (b) dynamical structure factors for various magnetizations 
m. Poles with a sizable weight represent stable magnons with a long lifetime whereas small weights indicate that spin-wave energies lie within 
a continuum. Magnons thus have a short lifetime and are susceptible to spontaneous decay. 



Away from the special (tt, tt) -point, our results are much 
clearer and show evidence for the decay of spin waves in multi 
spin- wave continua. Inspecting the longitudinal structure fac- 
tor at (37r/4, 37r/4) for instance, or equivalently, the trans- 
verse spin correlations at (7r/4,7r/4), we see that the spin- 
wave band fades away at m ^ 0.3 and enters a continuum, 
which on these finite clusters manifests itself as an area with 
densely distributed poles having relatively small intensity. For 



polarizations close to the ferromagnetic regime, the number 
of poles decreases rapidly. In the Hilbert space with only one 
flipped spin, the single pole coincides with the dashed classi- 
cal line because in the limit m 1/2, spin- wave results for 
the ferromagnet are exact. Compared to the situation at (tt, tt), 
which seems to be unique, our results at (37r/4, 37r/4) clearly 
show that the initially well defined magnon branch dives into 
a continuum and hence acquires a finite lifetime. 
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Looking at other points in the Brillouin zone, one clearly 
recognizes similar longitudinal instabilities at (7r/2,7r/2), 
(7r,0), and (7r/2,0), or equivalently at (7r,7r) shifted points 
in the transverse structure factor. The critical field at which 
magnon decay sets in depends on the momenta. To identify 
the regions of instability, it is useful to plot the weights of the 
first poles in the structure factors as a function of the magne- 
tization. Since different clusters have incompatible magneti- 
zation steps, it is necessary to use an interpolation to collect 
data from all clusters at a given magnetization. Figs. 15 show 
the weights of the lowest-lying pole in the longitudinal and 
the transverse structure factors for various magnetizations m. 
To eliminate spurious effects arising from very small Hilbert 
spaces at high polarization, we have only included data from 
sectors with m < mmax- We would like to emphasize that 
this way of illustrating the instabilities with respect to sponta- 
neous decay does not contain any information about the life- 
time of the magnons, but merely indicates the availability of 
phase space for decay. It is also important to note that in 
contrast to Fig. 13, representing the full structure factors, in 
Fig. 15, we only indicate the weights of the poles with low- 
est energy. This plot nicely illustrates the spreading out of the 
decay region, starting around (tt, tt) in the longitudinal and 
(0, 0) in the transverse structure factor, respectively. Note that 
this evolution is opposite to the instability evolution advocated 
in Ref. 20, where the decay into two spin waves sets in in the 
vicinity of the gapless point. We believe that this behavior 
is difficult to study, given the infrared cutoff imposed by the 
small size limitations of exact diagonalizations. On the other 
hand, our results clearly demonstrate the importance of multi- 
spin- wave continua for the decay processes, and those are first 
visible in the opposite comer of the Brillouin zone. 

Based on these results, we can say that high energy 
magnons around (tt, tt) become unstable at a critical polar- 
ization rric ~ 0.3. For higher magnetizations, the region of 
instability quickly increases. At m = 0.4, we expect the spin 
waves in more than half the Brillouin zone to have a finite 



lifetime. 



V. CONCLUSION 

We have studied the antiferromagnetic spin-^ Heisenberg 
model on the square lattice in a magnetic field. For the first 
time, we have calculated the field dependence of the mi- 
croscopic parameters of the low-energy long-wavelength a 
model description, namely the spin- wave velocity c, the spin- 
stiffness ps, and the transverse magnetic susceptibility x±- 
The motivation behind this work being the availability of com- 
pounds well described by this prototypical model, the main 
part of this paper has been devoted to studying the field de- 
pendence of the dynamical spin structure factors, directly pro- 
portional to the neutron scattering cross section. On the one 
hand, we hope that our comprehensive presentation of these 
dynamical quantities encourages future theoretical and espe- 
cially experimental work to examine particular aspects in even 
more detail. On the other hand, we are confident that our re- 
sults obtained by means of extensive exact diagonalizations 
of finite clusters with up to 64 sites, for which no approxima- 
tion or simplification has been made, serve as a benchmark for 
other methods. 
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